home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 1 / Precision Software Applications Silver Collection Volume One (PSM) (1993).iso / windows / games / xlmath21.arj / SOURCE.ZIP / XLMATH.C < prev    next >
C/C++ Source or Header  |  1993-04-10  |  21KB  |  851 lines
  1. /*----------------------------------------------------------------------*\
  2.  
  3.  XLMATH    
  4.  ========
  5.  Copyright    1992 by Roy Kari
  6.  
  7.  Author:    Roy Kari
  8.              Dept. of Chemistry
  9.             Laurentian University
  10.             Sudbury, Ont.
  11.             Canada        P3E 2C6
  12.             (705) 675-1151
  13.             Internet: "ROY@NICKEL.LAURENTIAN.CA"
  14.  
  15.  Revision history:
  16.  
  17.     1.0    September, 1992
  18.         initial version
  19.  
  20.     2.0    January, 1993
  21.         added dialog boxes
  22.  
  23.     2.1    April, 1993
  24.         bug fixes
  25.  
  26.  
  27.     This module contains routines written to interface normal C
  28.     routines to an Excel custom function format. All custom
  29.     functions are assumed to pass a floating point array to the custom
  30.     function, and the custom function returns an XLOPER, usually an
  31.     array. If you wish to use a dialog box, then the dialog routine
  32.     is written to interface the custom function to a dialog box.
  33.     The dialog interface must get the input from the Excel sheet,
  34.     change it to a floating point array, and then call the custom function.
  35.     On return from the custom function, the dialog interface must
  36.     display the return XLOPER in the sheet, and free all allocated
  37.     memory.
  38.  
  39.     2nd thoughts;
  40.     If I had this to over again, I would make the custom functions
  41.     accept their input in the form of XLOPER's rather than a 
  42.     floating point array. This would eliminate the need to convert.
  43.  
  44. \*----------------------------------------------------------------------*/
  45.  
  46.  
  47. /* --------------------------< Include files >--------------------- */
  48.  
  49. #include <windows.h>
  50. #include <xlcall.h>
  51.  
  52. #include <float.h>
  53. #include <math.h>
  54. #include <stdlib.h>
  55.  
  56. #include <framewrk.h>
  57. #include "xlmath.h"
  58. #include "xlutil.h"
  59. #include "xlmcurve.h"
  60. #include "xlmcfit.h"
  61.  
  62. REALTYPE predict (int k, PARVEC p, REALTYPE xk);
  63.  
  64. /********************************************************************
  65.  Diagonalize()
  66.  =============
  67.  Diagonalize a real symmetric matrix. Eigenvalues are returned
  68.  *******************************************************************/
  69. LPXLOPER PASCAL FAR __export Diagonalize(LPFP lpHmat)
  70. {
  71.     LPMULTI lpMulti = NULL;    // return XLOPER's
  72.     LPLPREAL lplpHmat;        // pointers for EXCEL allocated array
  73.     LPREAL lpEivec;            // eigenvector array
  74.     LPLPREAL lplpEivec;        // pointers for above
  75.     LPREAL lpEival;            // eigenvalue array
  76.  
  77.     WORD i,j,wRowIndex;        // index
  78.     int nIter;                // no iterations in _Diagonalize()
  79.     WORD wRows, wCols;        // saves typing
  80.     BOOL Success = TRUE;
  81.  
  82.     wRows = lpHmat->wRows;
  83.     wCols = lpHmat->wCols;
  84.  
  85.     // error if not square array
  86.     if (wRows != wCols)
  87.     {
  88.         ErrorHandler(XLM_NOT_SQUARE);
  89.         return (glpxlError);
  90.     }
  91.     
  92.     //     allocate a Global 2-D array for eigenvectors. This array allocation
  93.     //    demonstrates the use of 2-D array indexing M[i][j]. It is 
  94.     //    required in this case because the working routine _Diagonalize() was
  95.     //    written in this manner and it was too much work to convert
  96.     //    it to anything else.
  97.     if ((lpEivec = (LPREAL)GetMem( sizeof(double)*wRows*wCols)) == NULL)
  98.     {
  99.         return (glpxlError);
  100.  
  101.     }
  102.     if ((lplpEivec = InitPointers(lpEivec, wRows, wCols)) == NULL)
  103.     {
  104.         Success = FALSE;
  105.         goto Error3;        // free lpEivec
  106.     }
  107.  
  108.     // create pointers for Hmat; array allocated in Excel
  109.     if ((lplpHmat = InitPointers(&lpHmat->Data[0], wRows, wRows))
  110.             == NULL)
  111.  
  112.     {
  113.         Success = FALSE;
  114.         goto Error2;        // free lpEivec & lplpEivec
  115.     }
  116.  
  117.     // allocate a Global 1-D array for eigenvalues
  118.     if ((lpEival = (LPREAL)GetMem( sizeof(double)*wRows )) == NULL)
  119.     {
  120.         Success = FALSE;
  121.         goto Error1;        // free lpEivec, lplpEivec & lplpHmat
  122.     }
  123.  
  124.     // diagonalize the matrix
  125.     nIter = _Diagonalize(lplpHmat, lplpEivec, lpEival, wRows);
  126.  
  127.  
  128.     //
  129.     // return vales to EXCEL
  130.     //
  131.     
  132.     // init MULTI 1 row larger to store eigenvalues in last row
  133.     if ((lpMulti =InitMulti(wRows+1, wRows, xltypeNum)) == NULL)
  134.     {
  135.         Success = FALSE;
  136.     }
  137.     else
  138.     {
  139.         // copy the eigenvectors
  140.         for (i = 0; i < wRows; i++)
  141.         {
  142.             wRowIndex = i*wCols;
  143.             for (j=0; j < wCols; j++)
  144.             {
  145.                 lpMulti->Data[wRowIndex+j].val.num = lplpEivec[i][j];
  146.             }
  147.         }
  148.  
  149.         // copy the eigenvalues
  150.         wRowIndex = wRows*wCols;
  151.         for (j = 0; j < wCols; j++)
  152.         {
  153.             lpMulti->Data[wRowIndex+j].val.num = lpEival[j];
  154.         }
  155.     }
  156.  
  157.     // free pointers and arrays
  158.     FreeMem((LPVOID)lpEival);
  159. Error1:
  160.     MemFreePtr((LPVOID)lplpHmat);
  161. Error2:
  162.     MemFreePtr((LPVOID)lplpEivec);
  163. Error3:
  164.     FreeMem((LPVOID)lpEivec);
  165.     if (!Success)
  166.     {
  167.         return (glpxlError);
  168.     }
  169.  
  170.     // return XLOPER to EXCEL
  171.     return ( (LPXLOPER)lpMulti);
  172. }
  173.  
  174. /********************************************************************
  175.  MODensity()
  176.  =============
  177.  Compute the MO density as C(T)xOxC
  178.  *******************************************************************/
  179. LPXLOPER FAR PASCAL __export MODensity(LPFP lpEivec, LPFP lpOcc)
  180. {
  181.     LPMULTI lpDen = NULL;
  182.     UINT i, j, k;
  183.     UINT i_Idx, j_Idx, ij_Idx, NumVec;
  184.     BOOL bError = FALSE;
  185.  
  186.     NumVec = lpEivec->wRows;
  187.  
  188.     if (lpOcc->wRows == 1)
  189.     {
  190.         // Occ is row vector (1xN)
  191.         if (lpOcc->wCols != NumVec)
  192.         {
  193.             bError = TRUE;
  194.         }
  195.     }
  196.     else
  197.     {
  198.         // Occ is a column vector (Nx1)
  199.         if ( (lpOcc->wRows != NumVec) || (lpOcc->wCols != 1) )
  200.         {
  201.             bError = TRUE;
  202.         }
  203.     }
  204.     if (bError)
  205.     {
  206.         ErrorHandler(XLM_MODENSITY);
  207.         return (glpxlError);
  208.     }
  209.  
  210.     if ((lpDen =InitMulti(NumVec, NumVec, xltypeNum)) == NULL)
  211.     {
  212.         return (glpxlError);
  213.     }
  214.  
  215.     i_Idx = 0;
  216.  
  217.     for (i = 0; i < NumVec; i++)
  218.     {
  219.         j_Idx = 0;
  220.         for (j = 0; j < NumVec; j++)
  221.         {
  222.             ij_Idx = i_Idx+j;
  223.             lpDen->Data[ij_Idx].val.num = 0.0;
  224.             for (k = 0; k < NumVec; k++)
  225.             {
  226.                 // den[i][j] = Eivec[i][k]*Eivec[j][k]*Occ[k]
  227.                 lpDen->Data[ij_Idx].val.num += lpEivec->Data[i_Idx+k] *
  228.                                                 lpEivec->Data[j_Idx+k] *
  229.                                                 lpOcc->Data[k];
  230.             }
  231.             j_Idx += NumVec;
  232.         }
  233.         i_Idx += NumVec;
  234.     }
  235.     return ((LPXLOPER)lpDen);
  236. }
  237.  
  238.  
  239. /********************************************************************
  240.  PolyCurveFit()
  241.  =============
  242.  This function is a generalized least squares curve fitting function.
  243.  It fits a polynomial with linear coefficients to a dependent -
  244.  independent variable set of data
  245.  *******************************************************************/
  246. LPXLOPER PASCAL FAR _export PolyCurveFit(LPFP lpIndVar, LPFP lpDepVar,
  247.                 unsigned int Order)
  248. {
  249. #define    nCols    3    // placed in a NumObs x nCols array
  250.  
  251.     LPMULTI lpMulti;
  252.     LPREAL     lpYest, lpResid, lpPolycoef, lpCoefsig;
  253.     static double    See, Rsqrval, Rval, dferror;
  254.     static WORD wNumObs;
  255.     WORD wRows = lpIndVar->wRows;
  256.     WORD wCols = lpIndVar->wCols;
  257.     WORD wMinRows = 2*(Order+1) + 4;    // needed for return values
  258.     BOOL Success = TRUE;
  259.     BOOL bColumnVector = TRUE;
  260.     WORD i, wRow1, wRow2, wRowIndex;
  261.  
  262.     wNumObs = wRows;
  263.     if (wCols > wRows)
  264.     {
  265.         // change to row vectors
  266.         bColumnVector = FALSE;
  267.         wNumObs = wCols;
  268.     }
  269.  
  270.     // check input dimensions so both same
  271.     if     (     (lpDepVar->wRows != lpIndVar->wRows)     ||
  272.             (lpDepVar->wCols != lpIndVar->wCols)    ||
  273.             (Order >= wNumObs - 1)
  274.         )
  275.     {
  276.         ErrorHandler(XLM_POLY);
  277.         return(glpxlError);
  278.     }
  279.  
  280.     if (wMinRows < wNumObs)
  281.         wMinRows = wNumObs;
  282.  
  283.  
  284.     // memory for estimated Y values
  285.     lpYest = (LPREAL)GetMem( wNumObs*sizeof(double));
  286.     if (lpYest == NULL)
  287.     {
  288.         Success = FALSE;
  289.         goto Error4;
  290.     }
  291.     // memory for residuals
  292.     if ((lpResid = (LPREAL)GetMem( wNumObs*sizeof(double))) == NULL)
  293.     {
  294.         Success = FALSE;
  295.         goto Error3;
  296.     }
  297.  
  298.     //memory for coefficients of fitted polynomial
  299.     if ((lpPolycoef = (LPREAL)GetMem( (Order+1)*sizeof(double))) == NULL)
  300.     {
  301.         Success = FALSE;
  302.         goto Error2;
  303.     }
  304.  
  305.     // memory for standard errors of coefficient estimates
  306.     if ((lpCoefsig = (LPREAL)GetMem( (Order+1)*sizeof(double))) == NULL)
  307.     {
  308.         Success = FALSE;
  309.         goto Error1;
  310.     }
  311.  
  312.     Success  = _PolyCurveFit((LPREAL)&lpIndVar->Data[0],
  313.         (LPREAL)&lpDepVar->Data[0], (int)wNumObs, (int)Order, lpPolycoef,
  314.         lpYest, lpResid, (NPREAL)&See, lpCoefsig, &Rsqrval, &Rval, &dferror);
  315.     if (!Success)
  316.         goto Error0;
  317.  
  318.     // init XLOPER
  319.     if (bColumnVector)
  320.         lpMulti = InitMulti(wMinRows, nCols, xltypeNum);
  321.     else
  322.         lpMulti = InitMulti(nCols, wMinRows, xltypeNum);
  323.     if (lpMulti==NULL)
  324.     {
  325.         Success = FALSE;
  326.         goto Error0;
  327.     }
  328.  
  329.     // stuff return data into Multi
  330.     if (bColumnVector)
  331.     {
  332.         // return as column vectors (NumObs x 3)
  333.         wRowIndex = 0;
  334.         wRow2 = (Order+1)*nCols+2;
  335.            for (i=0; i < wNumObs;  i++)
  336.            {
  337.                lpMulti->Data[wRowIndex].val.num = lpYest[i];
  338.                lpMulti->Data[wRowIndex+1].val.num = lpResid[i];
  339.                if ( i <= Order )
  340.                {
  341.                    lpMulti->Data[wRowIndex+2].val.num = lpPolycoef[i];
  342.                    lpMulti->Data[wRowIndex+wRow2].val.num = lpCoefsig[i];
  343.                }
  344.             wRowIndex += nCols;
  345.            }
  346.            // remaining variables stuck in array directly below polycoef
  347.            wRowIndex = 2*(Order+1);
  348.         lpMulti->Data[(wRowIndex  )*nCols + 2].val.num = See;
  349.         lpMulti->Data[(wRowIndex+1)*nCols + 2].val.num = Rsqrval;
  350.         lpMulti->Data[(wRowIndex+2)*nCols + 2].val.num = Rval;
  351.         lpMulti->Data[(wRowIndex+3)*nCols + 2].val.num = dferror;
  352.  
  353.     }
  354.     else
  355.     {
  356.         // return as row vectors ( 3 x NumObs)
  357.         wRow1 = wNumObs;    // index to row 1
  358.         wRow2 = wNumObs*2;    // index to row 2
  359.         for (i=0; i<wNumObs; i++)
  360.         {
  361.             lpMulti->Data[i].val.num = lpYest[i];
  362.             lpMulti->Data[wRow1].val.num = lpResid[i];
  363.             ++wRow1;
  364.             if ( i<= Order)
  365.             {
  366.                 lpMulti->Data[wRow2].val.num = lpPolycoef[i];
  367.                 lpMulti->Data[wRow2+Order+1].val.num = lpCoefsig[i];
  368.                 ++wRow2;
  369.             }
  370.  
  371.         }
  372.         // remaining variables
  373.         wRow2 = wNumObs*2 + 2*(Order+1);
  374.         lpMulti->Data[wRow2].val.num = See;
  375.         lpMulti->Data[wRow2+1].val.num = Rsqrval;
  376.         lpMulti->Data[wRow2+2].val.num = Rval;
  377.         lpMulti->Data[wRow2+3].val.num = dferror;
  378.  
  379.     }
  380. Error0:
  381.     FreeMem((LPVOID)lpCoefsig);
  382. Error1:
  383.     FreeMem((LPVOID)lpPolycoef);
  384. Error2:
  385.     FreeMem((LPVOID)lpResid);
  386. Error3:
  387.     FreeMem((LPVOID)lpYest);
  388. Error4:
  389.     if (!Success)
  390.     {
  391.         return (glpxlError);
  392.     }
  393.     return ((LPXLOPER)lpMulti);
  394. }
  395. #undef nCols
  396.  
  397. /********************************************************************
  398.  CubicSplines()
  399.  =============
  400.  This function will fit a set of cubic spline polynomial equations
  401.  to a discrete set of data.
  402.  *******************************************************************/
  403. LPXLOPER PASCAL FAR _export CubicSplines(LPFP lpIndVar, LPFP lpDepVar)
  404. {
  405. #define    nCols    4    // placed in a wNumDat-1 x 4 array
  406.  
  407.     LPMULTI lpMulti;
  408.     WORD i, j, wIndex;
  409.     WORD wNumDat = lpIndVar->wRows;    // assume column vector
  410.     BOOL Success = TRUE;
  411.     LPREAL lpCoef;
  412.  
  413.     // check inputs are row or column vectors
  414.     if (    (lpIndVar->wCols > 1 && lpIndVar->wRows > 1)    ||
  415.             (lpDepVar->wCols > 1 && lpDepVar->wRows > 1) )
  416.         {
  417.             ErrorHandler(XLM_CUBIC);
  418.             return (glpxlError);
  419.         }
  420.  
  421.     // check if row or column input
  422.     if (lpIndVar->wRows < lpIndVar->wCols)
  423.     {
  424.         // input is a row vector
  425.         wNumDat = lpIndVar->wCols;
  426.     }
  427.  
  428.     if ((lpCoef = (LPREAL)GetMem( wNumDat*4*sizeof(realtype))) == NULL)
  429.         return (glpxlError);
  430.  
  431.     Success  = _CubicSplines((LPREAL)&lpIndVar->Data[0],
  432.                             (LPREAL)&lpDepVar->Data[0],
  433.                             (int)(wNumDat-1),
  434.                             (LPREAL)lpCoef);
  435.     if (!Success)
  436.     {
  437.         FreeMem((LPVOID)lpCoef);
  438.         ErrorHandler(XLM_CUBIC);
  439.         return (glpxlError);
  440.     }
  441.  
  442.     // always return coefs in N-1 x 3 matrix
  443.     if ((lpMulti = InitMulti(wNumDat-1, nCols, xltypeNum)) == NULL)
  444.     {
  445.         FreeMem((LPVOID)lpCoef);
  446.         return (glpxlError);
  447.     }
  448.     for (i=0; i < wNumDat-1; i++)
  449.     {
  450.         wIndex = i*nCols;
  451.         for (j = 0; j < nCols; j++)
  452.         {
  453.             lpMulti->Data[wIndex + j].val.num = lpCoef[wIndex + j];
  454.         }
  455.     }
  456.     FreeMem((LPVOID)lpCoef);
  457.     return ((LPXLOPER)lpMulti);
  458. }
  459. #undef nCols
  460.  
  461. /********************************************************************
  462.  CalcSpline()
  463.  =============
  464.  This function calculate the cubic spline interpolation of an
  465.  y-value given an x-value and the cubic spline coefficients.
  466.  
  467.  *******************************************************************/
  468. LPXLOPER PASCAL FAR __export CalcSpline(LPFP lpXorig, LPFP lpCoef,
  469.                                             LPFP lpXcalc)
  470. {
  471.     LPMULTI lpMulti;
  472.     WORD wRows;
  473.     WORD wCols;
  474.     WORD wNumCalcX;
  475.     WORD wNumOrigX;
  476.     realtype Xmin, Xmax, Xcurrent;
  477.     WORD i;
  478.  
  479.     // check orientation of input vectors
  480.     if (lpXorig->wRows > lpXorig->wCols)
  481.     {
  482.         // column vector
  483.         wNumOrigX = lpXorig->wRows;
  484.         wRows = wNumCalcX = lpXcalc->wRows;
  485.         wCols = 1;
  486.     }
  487.     else
  488.     {
  489.         // row vector
  490.         wNumOrigX = lpXorig->wCols;
  491.         wCols = wNumCalcX = lpXcalc->wCols;
  492.         wRows = 1;
  493.     }
  494.  
  495.     // check limits
  496.     if     (     (lpCoef->wRows != wNumOrigX - 1) ||
  497.              (lpCoef->wCols != 4)
  498.         )
  499.     {
  500.         ErrorHandler(XLM_CALCSPLINE);
  501.         return glpxlError;
  502.     }
  503.  
  504.     // allocate memory for returned xloper
  505.     if ( (lpMulti = InitMulti(wRows, wCols, xltypeNum)) == NULL)
  506.         return (glpxlError);
  507.  
  508.     // do the calculations
  509.     wNumOrigX = wNumOrigX - 1;;
  510.     Xmin = lpXorig->Data[0];
  511.     Xmax = lpXorig->Data[wNumOrigX];
  512.  
  513.     for (i=0; i < wNumCalcX; i++)
  514.     {
  515.         Xcurrent = lpXcalc->Data[i];
  516.         if ( (Xcurrent < Xmin) || (Xcurrent > Xmax))
  517.         {
  518.             lpMulti->Data[i].xltype = xltypeErr;
  519.             lpMulti->Data[i].val.err = xlerrValue;
  520.         }
  521.         else
  522.         {
  523.         _CalcSpline(    (LPREAL)&lpXorig->Data[0],
  524.                         (LPREAL)&lpCoef->Data[0],
  525.                         (int)(wNumOrigX),
  526.                         (LPREAL)&Xcurrent,
  527.                         (LPREAL) &lpMulti->Data[i].val.num);
  528.         }
  529.     }
  530.     return ((LPXLOPER)lpMulti);
  531. }
  532. /********************************************************************
  533.  SmoothSG()
  534.  =============
  535.  This function uses the Savitsky - Golay algorithm to reduce the
  536.  noise in a sampled data set.
  537.   *******************************************************************/
  538. LPXLOPER PASCAL FAR _export SmoothSG(LPFP lpData, unsigned int wSmoothNum,
  539.     unsigned int wDerivNum)
  540. {
  541.     LPMULTI lpMulti;
  542.     WORD wRows = lpData->wRows;
  543.     WORD wCols = lpData->wCols;
  544.     WORD wNumDat;
  545.     LPREAL lpCoef;
  546.     WORD i;
  547.  
  548.     // check limits
  549.     if ( (wSmoothNum < 1) || (wSmoothNum > 5) || (wDerivNum > 2) )
  550.     {
  551.         ErrorHandler(XLM_SMOOTHSG);
  552.         return glpxlError;
  553.     }
  554.  
  555.     // set wNumDat & check if row or column vector
  556.     wNumDat = wRows > wCols? wRows : wCols;
  557.     if (wCols > 1 && wRows > 1)
  558.         {
  559.             ErrorHandler(XLM_SMOOTHSG);
  560.             return (glpxlError);
  561.         }
  562.  
  563.     if ((lpCoef = (LPREAL)GetMem(wNumDat*sizeof(double))) == NULL)
  564.         return glpxlError;
  565.  
  566.     // call smoothing routine
  567.     _DataSmoothSg((LPREAL)&lpData->Data[0],
  568.                     (int)wNumDat,
  569.                     (int)wSmoothNum,
  570.                     (int)wDerivNum,
  571.                     (LPREAL)lpCoef);
  572.  
  573.     // copy results & return as row or column vector in original
  574.     if ((lpMulti = InitMulti(wRows, wCols, xltypeNum)) == NULL)
  575.     {
  576.         FreeMem((LPVOID)lpCoef);
  577.         return glpxlError;
  578.     }
  579.     else
  580.     {
  581.         for (i=0; i<wNumDat; i++)
  582.         {
  583.             lpMulti->Data[i].val.num = lpCoef[i];
  584.         }
  585.     }
  586.     FreeMem((LPVOID)lpCoef);
  587.  
  588.     return ((LPXLOPER)lpMulti);
  589. }
  590. /********************************************************************
  591.  SmoothWeights()
  592.  =============
  593.  This function smooths data via a weighted average
  594.  *******************************************************************/
  595. LPXLOPER PASCAL FAR _export SmoothWeights(LPFP lpData,
  596.             LPFP lpWeights)
  597. {
  598.     LPMULTI lpMulti;
  599.     WORD wRows = lpData->wRows;
  600.     WORD wCols = lpData->wCols;
  601.     WORD wNumDat;
  602.     WORD wSmoothNum;
  603.     LPREAL lpCoef;
  604.     realtype Divisor = 0.0;
  605.     WORD i;
  606.  
  607.     // set wNumDat & check if row or column vector
  608.     wNumDat = wRows > wCols? wRows : wCols;
  609.     if (wCols > 1 && wRows > 1)
  610.         {
  611.             ErrorHandler(XLM_SMOOTHSG);
  612.             return (glpxlError);
  613.         }
  614.  
  615.     // take from column or row vector
  616.     wSmoothNum = (lpWeights->wRows > lpWeights->wCols) ? lpWeights->wRows :
  617.                     lpWeights->wCols;
  618.  
  619.     // get normalization factor for weights
  620.     for (i=0; i<wSmoothNum; i++)
  621.         Divisor += lpWeights->Data[i];
  622.  
  623.     if ((lpCoef = (LPREAL)GetMem(sizeof(double)*wNumDat)) == NULL)
  624.         return glpxlError;
  625.  
  626.     _DataSmoothWeights((LPREAL)&lpData->Data[0],
  627.                     (int)wNumDat,
  628.                     (int)wSmoothNum,
  629.                     (LPREAL)&lpWeights->Data[0],
  630.                     (LPREAL) &Divisor,
  631.                     (LPREAL)lpCoef);
  632.  
  633.     // return as column or row vector as per original
  634.     if ((lpMulti = InitMulti(wRows, wCols, xltypeNum)) == NULL)
  635.     {
  636.         FreeMem((LPVOID)lpCoef);
  637.         return glpxlError;
  638.     }
  639.     else
  640.     {
  641.         for (i=0; i<wNumDat; i++)
  642.         {
  643.             lpMulti->Data[i].val.num = lpCoef[i];
  644.         }
  645.     }
  646.     FreeMem((LPVOID)lpCoef);
  647.     return ((LPXLOPER)lpMulti);
  648. }
  649. /********************************************************************
  650.  CustomFit()
  651.  =============
  652.  This function is the interface function to _CustomFit()
  653.  *******************************************************************/
  654. LPXLOPER PASCAL FAR __export CustomFit(LPFP lpData, LPFP lpParms)
  655. {
  656.     extern LPMULTI lpMultiCfit;
  657.     LPMULTI lpMulti;
  658.     WORD wRows = lpData->wRows;
  659.     WORD wCols = lpData->wCols;
  660.     static WORD wNumDat;
  661.     static WORD wNumParm;
  662.     LPREAL lpY, lpX, lpP, lpUpper, lpLower;
  663.     LPREAL lpPopt, lpSdev, lpD, lpEig, lpFit;
  664.     WORD i, wRowIndex;
  665.     BOOL success = FALSE;
  666.  
  667.     // set wNumDat & check if row or column vector
  668.     if (wRows > wCols)
  669.     {
  670.         // in column vector -- OK
  671.         wNumDat = wRows;
  672.         if (wCols !=2 || wRows <3)
  673.         {
  674.             // need 2 cols & more than 2 rows
  675.             ErrorHandler(XLM_CF1);
  676.             return(glpxlError);
  677.  
  678.         }
  679.     }
  680.     else
  681.     {
  682.         // in row vector & NOT allowed
  683.         ErrorHandler(XLM_CF2);
  684.         return(glpxlError);
  685.     }
  686.  
  687.     // Parms must be row vector
  688.     wNumParm =     lpParms->wCols;
  689.     if (lpParms->wRows != 3)
  690.     {
  691.         // must have 3 rows
  692.         ErrorHandler(XLM_CF3);
  693.         return(glpxlError);
  694.         
  695.     }
  696.  
  697.     // must have more data points than parms
  698.     if (wNumParm > wNumDat)
  699.     {
  700.         ErrorHandler(XLM_CF8);
  701.         return(glpxlError);
  702.     }
  703.     // allocate memory
  704.     if ((lpY = (LPREAL)GetMem(sizeof(double)*wNumDat)) == NULL)
  705.         return glpxlError;
  706.     if ((lpX = (LPREAL)GetMem(sizeof(double)*wNumDat)) == NULL)
  707.     {    
  708.         FreeMem((LPVOID)lpY);
  709.         return glpxlError;
  710.     }
  711.     if ((lpP = (LPREAL)GetMem(sizeof(double)*wNumParm)) == NULL)
  712.     {
  713.         FreeMem((LPVOID)lpY);
  714.         FreeMem((LPVOID)lpX);
  715.         return glpxlError;
  716.     }
  717.     if ((lpUpper = (LPREAL)GetMem(sizeof(double)*wNumParm)) == NULL)
  718.     {
  719.         FreeMem((LPVOID)lpY);
  720.         FreeMem((LPVOID)lpX);
  721.         FreeMem((LPVOID)lpP);
  722.         return glpxlError;
  723.     }
  724.     if ((lpLower = (LPREAL)GetMem(sizeof(double)*wNumParm)) == NULL)
  725.     {
  726.         FreeMem((LPVOID)lpY);
  727.         FreeMem((LPVOID)lpX);
  728.         FreeMem((LPVOID)lpP);
  729.         FreeMem((LPVOID)lpUpper);
  730.         return glpxlError;
  731.     }
  732.     
  733.     // extract the data
  734.     wRowIndex = 0;
  735.     for (i=0; i<2*wNumDat; i+=2)
  736.     {
  737.         lpX[wRowIndex] = lpData->Data[i];
  738.         lpY[wRowIndex] = lpData->Data[i+1];
  739.         ++wRowIndex;
  740.     }
  741.  
  742.     // extract the parms
  743.     for (i=0; i<wNumParm; i++)
  744.     {
  745.         lpP[i]         = lpParms->Data[i];
  746.         lpUpper[i]     = lpParms->Data[i+wNumParm];
  747.         lpLower[i]     = lpParms->Data[i+2*wNumParm];
  748.     }
  749.  
  750.     // allocate fitting variable memory
  751.     if ((lpPopt = (LPREAL)GetMem(sizeof(double)*wNumParm)) == NULL)
  752.     {    
  753.         FreeMem((LPVOID)lpPopt);
  754.         success = FALSE;
  755.         goto abort;
  756.     }
  757.     if ((lpSdev = (LPREAL)GetMem(sizeof(double)*wNumParm)) == NULL)
  758.     {    
  759.         FreeMem((LPVOID)lpPopt);
  760.         success = FALSE;
  761.         goto abort;
  762.     }
  763.     if ((lpD = (LPREAL)GetMem(sizeof(double)*wNumParm)) == NULL)
  764.     {    
  765.         FreeMem((LPVOID)lpPopt);
  766.         FreeMem((LPVOID)lpSdev);
  767.         success = FALSE;
  768.         goto abort;
  769.  
  770.     }
  771.     if ((lpEig = (LPREAL)GetMem(sizeof(double)*wNumParm)) == NULL)
  772.     {    
  773.         FreeMem((LPVOID)lpPopt);
  774.         FreeMem((LPVOID)lpSdev);
  775.         FreeMem((LPVOID)lpD);
  776.         success = FALSE;
  777.         goto abort;
  778.     }
  779.     if ((lpFit = (LPREAL)GetMem(sizeof(double)*wNumParm)) == NULL)
  780.     {    
  781.         FreeMem((LPVOID)lpPopt);
  782.         FreeMem((LPVOID)lpSdev);
  783.         FreeMem((LPVOID)lpD);
  784.         FreeMem((LPVOID)lpEig);
  785.         success = FALSE;
  786.         goto abort;
  787.     }
  788.  
  789.     // allocate a multi for calling macro
  790.     // strip xlbitDLLFree to prevent Excel from freeing 
  791.     if ((lpMultiCfit = InitMulti(1, wNumParm, xltypeNum)) == NULL)
  792.         goto abort1;
  793.     lpMultiCfit = (LPMULTI)StripxlbitDLLFree((LPXLOPER)lpMultiCfit);
  794.     
  795.     
  796.     success = _CustomFit((int)wNumDat, (int)wNumParm, 
  797.                             lpY, lpX, lpP, lpUpper, lpLower,
  798.                             lpPopt, lpSdev, lpD, lpEig, lpFit);
  799.     if (!success)
  800.     {
  801.         goto abort2;        
  802.     }
  803.     //
  804.     // do output via Multi
  805.     // 
  806.     if ((lpMulti = InitMulti(wNumDat, wNumParm+1, xltypeNum)) == NULL)
  807.     {
  808.         goto abort2;
  809.     }
  810.  
  811.     // calculate fitted values to first column
  812.     wRowIndex = 0;
  813.     for (i=0; i<wNumDat; i++)
  814.     {
  815.         lpMulti->Data[wRowIndex].val.num = predict(i, lpPopt, lpX[i]);
  816.         wRowIndex +=(wNumParm+1);
  817.     }
  818.  
  819.     // copy parameters, std. deviations & eigenvalues
  820.     wRowIndex = wNumParm+1;
  821.     for (i=0; i<wNumParm; i++)
  822.     {
  823.         lpMulti->Data[i+1].val.num                 = lpPopt[i];
  824.         lpMulti->Data[wRowIndex+i+1].val.num     = lpSdev[i];
  825.         lpMulti->Data[2*wRowIndex+i+1].val.num     = lpEig[i];
  826.         lpMulti->Data[3*wRowIndex+i+1].val.num    = lpD[i];
  827.     }
  828. abort2:
  829.     xlAutoFree((LPXLOPER)lpMultiCfit);
  830.     Excel(xlFree, 0, 1, (LPXLOPER)&xMacroRef);
  831. abort1:
  832.     FreeMem(lpFit);
  833.     FreeMem(lpEig);
  834.     FreeMem(lpD);
  835.     FreeMem(lpSdev);
  836.     FreeMem(lpPopt);
  837.  
  838. abort:    
  839.     FreeMem((LPVOID)lpY);
  840.     FreeMem((LPVOID)lpX);
  841.     FreeMem((LPVOID)lpP);
  842.     FreeMem((LPVOID)lpUpper);
  843.     FreeMem((LPVOID)lpLower);
  844.     if (success)
  845.         return ((LPXLOPER)lpMulti);
  846.     else
  847.     {
  848.         return(glpxlError);
  849.     }
  850. }    
  851.